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ABSTRACT 

We present the first global, three-dimensional simulations of solar/stellar convection that take into 
account the influence of magnetic flux emergence by means of the Babcock-Leighton (BL) mechanism. 
We have shown that the inclusion of a BL poloidal source term in a convection simulation can pro- 
mote cyclic activity in an otherwise steady dynamo. Some cycle properties are reminiscent of solar 
observations, such as the equatorward propagation of toroidal flux near the base of the convection 
zone. However, the cycle period in this young sun (rotating three times faster than the solar rate) is 
very short (~ 6 months) and it is unclear whether much longer cycles may be achieved within this 
modeling framework, given the high efficiency of field generation and transport by the convection. 
Even so, the incorporation of mean-field parameterizations in 3D convection simulations to account 
for elusive processes such as flux emergence may well prove useful in the future modeling of solar and 
stellar activity cycles. 

Subject headings: Sun:dynamo, Sun:interior, Stars: magnetic field, Convection, Magnetohydrodynam- 
ics (MHD) 



1. INTRODUCTION 

The emergence of magnetic flux through the solar 
photosphere regulates solar variability and powers space 
weather. It is clear that this flux originates in the solar 
interior and is produced by the solar dynamo. However, 
it is currently unclear what role flux emergence plays in 
establishing the 22-year solar activity cycle. Is it an es- 
sential ingredient or merely a superficial by-product of 
deeper-seated dynamics? 

One of the principle means by which flux emergence 
may act an an essential ingredient in the operation of the 
solar dynamo is t hrough the so-called Babcock-Le ighton 
(BL) mechanism (|Babcock| |196li; iLei ghtoni fl964) . The 
BL mechanism arises through the dynamics of flux emer- 
gence. As a buoyant flux tube rises through the con- 
vection zone, the Coriolis force induces a twist in the 
axis of the tube that is manifested upon emergence as 
a poleward displacement between the trailing and lead- 
ing edges of a bipolar active region. When the active 
region subsequently fragments and disperses due to sur- 
face convection and meridional flow, the redistribution of 
vertical flux induces a net electromotive force (emf) that 
converts mean toroidal field to poloidal field. Although 
doubts remain about the viability of the BL mechanism 
as the principle source of poloidal flux, its empirical foun- 
dation and robustness have made it an integral element 
in many recent m ean- field dynamo models of the solar cy- 
cle (r eviewed by iDikpati fc Gilmanl 120091 : iCharbonneaul 

[2oTm . 

Current simulations of global solar and stellar convec- 
tion do exhibit sustained dynamo action and, depending 
on the parame t er reg i me, can prod uce ma gnetic cycles 
dGhizaru et al.l l201Ct iRacine et all 120111: Brown et ahl 
1201 lh . Yet, these simulations do not have sufficient res- 
olution to r eliably captur e flux emergence or the BL 
mechanism. [Nelson et al.l f|201 lh recently reported the 
first convective dynamo simulation to exhibit the spon- 



taneous, self-consistent generation of buoyant magnetic 
flux structures generated by convectively-driven rota- 
tional shear. However, even these cannot realistically 
simulate the subtle, multi-scale dynamics of flux emer- 
gence and dispersal that underlies the BL mechanism. 
This would require (1) very low diffusion to form concen- 
trated flux structures, (2) very high resolution to capture 
the destabilization and coherent rise of those structures, 
and (3) a realistic depiction of surface convection, merid- 
ional flow, and radiative transfer in order to capture 
the emergence, coalescence, fragmentation , and disper- 
sal of bipolar active regions (jCheung et al.l[2010h . Each 
of these is a formidable computational challenge in its 
own right, stetching the limits of modern supercomput- 
ers. No global model is capable of unifying convective 
dynamos and the BL mechanism through direct numeri- 
cal simulation. 

Thus, the presence of magnetic cycles in global convec- 
tion simulations demonstrates that the BL mechanism is 
not a necessary ingredient for cyclic activity. How this 
may or may not apply to the solar dynamo is a complex, 
unresolved issue and we make no attempt at a compre- 
hensive discussion here. Our purpose is rather to inves- 
tigate how flux emergence may alter the behavior of a 
convective dynamo by means of the BL mechanism. The 
BL mechanism is modeled using a mean-field parame- 
terization intended to mimic dynamics that cannot be 
explicitly captured by the simulation itself for the rea- 
sons outlined above. Our approach is described in $2] 
and simulation results are presented in ^J3j along with 
interpretive discussion. 

Before proceeding, it is worth emphasizing up front 
that the simulations we consider here are of a solar- 
like star rotating three times more rapidly than the Sun 
(3f2©). This is done because it is a tidy numerical ex- 
periment; without BL forcing, this dynamo builds strong 
large-scale fields that do not undergo cycles. Other pa- 



2 



Miesch & Brown 




-50 



4 6 
time (years) 



10 



Fig. 1. — Progenitor simulation, Case D3, from [Brown et al.l i20Wt ). (a) Radial velocity in a Mollweide projection near the top of the 
computational domain (r = 0.953-R, t = 2 yr, color table saturation ± 70 m s — 1 , yellow/orange upflow, blue/black downflow). (b) Angular 
velocity Q/2ir, averaged over 1200 days (7.8-11.3 yrs, saturation 1140-1320 nHz, white/red fast, blue/black slow), (c) Mean toroidal 
magnetic field (-£?<;>) averaged over longitude and time (15 day average near t = 2 yrs, saturation ± 8 kG, red positive, blue negative), (d) 
(-B</>) in the mid convection zone (r = 0.84i?) versus latitude and time (colors as in c, saturation ± 4 kG). 



rameter regimes, including those at the solar rotation 
rate, will be considered in future work. 

2. MODEL 

The starting point for our investigation is the convec- 
tive dynamo simulation that we refer to as c ase D 3 and 
that is described in detail by IBrown et al.l (|2010D . We 
refer the reader to that paper for further information on 
the set up and results of the simulation as well as the 
ASH (Anelastic Spherical Harmonic) code that is used 
to solve the equations of magnetohydrodynamics (MHD) 
in a rotating spherical shell under the anelastic approxi- 
mation. Solar values are used for the luminosity and the 
background stratification but the rotation rate is a factor 
of three faster than the Sun, as mentioned above. The 
spatial resolution of all simulations reported in this pa- 
per is 96 x 256 x 512 (r, 9, <fi). The simulation domain 
spans from r\ ~ 0.718-R to r 2 = 0.966i?, where R is the 
solar radius. 

The salient feature of Case D3 that we are inter- 
ested in here is the presence of persistent toroidal field 
structures that we term magnetic wreaths. The con- 
vection exhibits an intricate, evolving small-scale struc- 
ture (Fig. [TJi) but produces a substantial differential ro- 
tation (Fig. [lj) that in turn promotes the generation 
of prominent magnetic wreaths (Fig. UJp). The mean 
toroidal field is approximately symmetric about the equa- 
tor, with one wreath in each hemisphere of opposite po- 
larity. Although continual buffeting by convective mo- 
tions induces non-axisymmetric and temporal fluctua- 
tions, the mean field is remarkably persistent, retain- 
ing its essential structure indefinitely (Fig. Q]i). After 
they are established, the wreaths persist for at least 60 
years (t he duration of the simulation) with no sign of 
abating ([Brown et a l. 2010). This time interval is much 
longer than the rotation period (9.3 days), the convec- 
tive turnover time scale (about 20 days) and the ohmic 
diffusion time (about 3.6 years). 

The simulations described in this paper are all 
restarted from the same iteration of Case D3, defined as 
time t = 0. This includes the unmodified D3 run shown 
in Figure [T] which was continued beyond the restart it- 
eration for comparison purposes. The absence of initial 
transients in Fig.[T]demonstrates that this simulation has 
reached a statistically steady state by the mutual refer- 
ence time t = 0. 

The anelastic MHD equations for the conservation of 
mass, momentum, and thermal energy are solved with 



no modification. The only difference betwee n the simula- 
tions presented here and that presented in IBrown et al.l 
(2010) is in the magnetic induction equation where we 
add an additional poloidal source term S(r, 9) as follows: 



dB „ 



vxB - ??VxB + S<f>) 



(1) 



The additional term is intended to mimic the generation 
of mean poloidal field by the BL me chanism . Follo wing 
the mean- field BL dynamo model of Rempel (20Q6|), we 
choose 

S(r,9)=af(r)g(9)B 4> (2) 
where /(r) and g{9) are radial and latitudinal profiles 

(r - r 2 ) 
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d 2 
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3V3 



sin 2 9 cos 9 



(3) 



(4) 



and is a measure of the mean toroidal flux at the base 
of the convection zone 



/ > 



h(r) (B$) dr 



(5) 



Brackets <> denote an average over longitude and h is 
an averaging kernal given by h(r) = ho (r — r\) (r& — r). 
The integration is confined to a region near the base of 
the convection zone, below r\, = 0.79i? and the normal- 
ization ho is defined such that J r h(r)dr = 1. 

Note that the radial profile in equation ([3]) is nonzero 
only near the top of the convection zone, for r > r 2 — d. 
We choose d = 20 Mm so the poloidal source operates 
above r — 0.937i?. Thus, the BL term is nonlocal in the 
sense that the poloidal source near the surface is pro- 
portional to the mean toroidal flux near the base of the 
conv ection zone. T his is typical of BL dynamo models 
(e.g. lRempelll2006l ). 

The amplitude of the BL term, a, includes an al- 
gebraic quenching of the form a — a (l + B' 2 / B 2 )^ 1 
where B q = 1 MG is the quenching field strength and 
B 2 = (1/2) B 2 sin 9d9. In practice the fields gener- 
ated are much less than B q so the quenching plays little 
role. 

The magnetic diffusivity rj in all simulations is the same 
as in Case D3, varying from 1.56-7.69 xlO 12 cm 2 s _1 
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Fig. 2. — Mean toroidal field as a function of latitude and time 
as in Fig. [TJZ (same radial level, color scale, saturation ± 4 kG) 
with (a) ao = 1 m s~\ (6, c) ao = 10 m s —1 , and (<f, e) ao = 
100 m s — 1 . Frame (c) represents a continuation of frame (b) at 
much later times. Frame (e) is a zoomed-in portion of frame (d) 
highlighting the magnetic cycles. 

from the bottom to the top of the convection zone, cx 
p -1 / 2 , where p is the background density. 

The objective of this paper is to vary the amplitude 
of the BL term ao in order to investigate how flux emer- 
gence may alter the nature of the dynamo. We emphasize 
again that the 3D MHD equations are unmodified apart 
from the S(r, 9) term in eq. ([T]) so setting ao — repro- 
duces Ca se D3 as shown in F igure Q] and as described at 
length in lBrown et al.l ([2010D . 

3. RESULTS 

Figure [2] shows "butterfly" diagrams (latitude-time 
plots of the mean toroidal field near the base of the con- 
vection zone) for a series of simulations with progres- 



sively increasing values of the BL forcing amplitude ao- 
For c«o = 1 m s _1 the BL term is too weak to sigifi- 
cantly influence the operation of the dynamo (Fig. Wp). 
The axisymmteric component of the poloidal field near 
the surface does increase but not enough to effect the 
maintenance of persistent wreaths in the lower convec- 
tion zone. 

For «o = 10 m s _1 the results are very different. The 
wreaths essentially vanish within a few years of simula- 
tion time (Fig. Wi>)- This can be attributed to the dif- 
ferential rotation operating on the mean poloidal field 
generated by the BL mechanism via the fi-effect. The 
sense of the poloidal field produced by the BL term is 
such that the f2-effect generates toroidal flux in the up- 
per convection zone that is of opposite polarity to the 
sense of the wreaths. Convection rapidly mixes this flux, 
bringing together opposite polarities that are annihilated 
through ohmic dissipation. The magnetic energy in the 
mean fields drops rapidly, decaying by a factor of 10 6 by 
t ~ 30 yrs. 

Yet, beyond about 30 years, the magnetic energy be- 
gins to rise again and the dynamo is reborn, saturating 
by about t ~ 75 yrs (Fig. [5b) with a very different struc- 
ture. Persistent toroidal wreaths are again present but 
they are symmetric about the equator, with two wreaths 
per hemisphere and a somewhat lower magnetic energy 
(about 27% less than in Case D3). 

Increasing ao by another order of magnitude produces 
prominent magnetic cycles (Fig. ®i) . Toroidal wreaths 
form at mid latitudes in each hemisphere and propagate 
toward the equator, reminiscent of the solar butterfy dia- 
gram (Fig. [5b). However, the cycle period is much shorter 
than in the Sun; about 6 months compared to 22 years. 
Nonlinear modulation of the cycle is evident, with waxing 
and waning amplitudes and transient, weaker substruc- 
ture in the butterfly diagram near the equator. 

It is clear from Fig. [5b that reversals in the two hemi- 
spheres are not precisely sychronized. A more careful 
analysis reveals that the phase difference shifts over time, 
suggesting that the two hemispheres are to some extent 
decoupled. For example, the northern hemisphere leads 
the southern hemisphere from t ~ 4.5 - 11.5 yrs but 
the reverse is true for t ~ 12 - 14.5 yrs. Similar shifts in 
the hemispheric phase difference (albeit less pronounced) 
have been reported in sunspot records; in particular, 
the southern hemisphere was apparently leading in cy- 
cles 18-20 while the no rthern was leading in cycles 21-23 
([Mcintosh et al.ll20l"l . 

The symmetry of the dynamo about the equator can 
be quantified by the parity V = (B 2 - B 2 a )/(B 2 S + B 2 a ), 
where B s and B a are the symmetric and antisymmtric 
components of (B^) respectively (sampled at r = 0.84i?). 
The parity of Case D3 and for ao = 1 m s _1 ranges be- 
tween -0.5 and -1 while that for ao = 10 m s^ 1 ranges 
between positive 0.5-1. By contrast, the cyclic dynamo 
(ao = 100 m s _1 ) shifts between positive and negative 
parity as time proceeds. It does not exhibit the high syn- 
chronization of the solar cycle which exhibits a persistent 
negative parity. 

As expected, the addition of a BL term has a sub- 
stantial influence on the amplitude and structure of the 
mean poloidal field. Without it, the poloidal field has 
a roughly octupolar structure, such that (B r ) is radially 
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Fig. 3. — Parity of the mean toroidal field as a function of time 
for the simulation shown in Fig.[2]i,e (c»o = 100 m s _1 ). 

outward near the north pole, inward near the core of the 
northern wreath, a nd antisymmetric about the equator 
(|Brown et al.l 12010). As noted above, the BL term gen- 
erates opposing poloidal flux at mid-latitudes near the 
surface, producing multi-polar structure and enhancing 
dissipation. This plays an essential role in establishing 
the cycles of Fig [2]i-e and is evident in movies of the 
mean-field evolution. The typical amplitude of the mean 
poloidal field near the surface for ao = 100 m s _1 , 1-2 
kG, is much larger than for ao = (~ 200 G in Case D3) 
but the overall (3D) magnetic energy is about a factor of 
two smaller. 

The transition from steady to cyclic dynamos occurs 
when the BL term S competes with the fluctuating emf 
(v' XB'), where v' = v - (v) and B' = B - (B). This in 
turn occurs when ao becomes comparable to the velocity 
of the convective motions, V c . The relevant scale for V c is 
the rms value of the meridional components of v', which 
is about 120 m s" 1 near the top of the convection zone 
in Case D3, decreasing with depth. 

A legitimate question is whether the cyclic dynamo in 
Fig. [5Ji,e is operating in an essentially mean-field, ax- 
isymmteric mode. In other words, if the BL term were 
to dominate the generation of mean poloidal field and 
the differential rotation were to dominate the generation 
of mean toroidal field via the Si-effect, then one might 
expect this 3D simulation to behave very similarly to an 
analogous, axisymmetric mean-field model. Under this 
mean-field scenario, coupling between the poloidal and 
toroidal source regions might occur through the nonlo- 
cality of the BL term, the mean meridional flow, and the 
turbulent diffusion rj. The primary role of the resolved 
convective motions would then be to maintain the mean 
flows. 

In order to assess whether this mean-field scenario is in- 
deed a valid interpretation of the cyclic activity shown in 
Fig. ERe, we have initiated another simulation in which 
we have artificially suppressed the fluctuating emf. More 
specifically, we have replaced the v X B term in equation 
(fT]) with only mean-field induction (v) X (B). This sim- 
ulation was restarted from that shown in Fig. [2}i,e at 
t = 8.25 yrs with the same parameters. The only differ- 
ence is the absence of the fluctuating emf. 

The dynamo mode changes dramatically, as demon- 
strated in Fig. [H The non-axisymmctric magnetic field 
is quickly dissipated by ohmic diffusion, with the corre- 
sponding energy decreasing by six orders of magnitude 
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Fig. 4. — Mean toroidal field versus latitude and time as in Fig. [2] 
(same radial level, color table, saturation ± 4kG) for a simulation 
with ao = 100 m s — x , along with an artificial suppression of the 
fluctuating emf (v' X B') (see text). Compare with Fig. [2]i,e. 

within two years (shear promotes more rapid dissipation 
than the nominal 3.6 year diffusion time scale). The total 
magnetic energy is about a factor of three larger than in 
the progenitor simulation of Fig. [2]i,e and is dominated 
by a strong, axisymmtric toroidal field that is symmetric 
about the equator (positive parity) and reverses cycli- 
cally with a period of about 1.3 yrs. The butterfly dia- 
gram in Fig. U suggests poleward propagation but closer 
inspection of the mean field evolution reveals a dynamo 
wave that propagates toward the rotation axis, with a 
cylindrical orientation for the wave front. This is consis- 
tent with the cylindrical nature of the differential rota- 
tion profile (Fig. [!}>) but is strikingly different from the 
progenitor simulation in Fig. [2Ji,e which exhibits virtu- 
ally no mean toroidal field at the equator even when the 
parity is positive. In short, cycles are longer, stronger, 
and less solar-like without a turbulent emf. 

Thus, the cyclic dynamo in Fig. [2ji,e is not operat- 
ing as a mean-field dynamo, or at least not in a naive 
sense. Convection contributes mean field generation and 
transport that plays an essential role in shaping the dy- 
namo even for ao as high as 100 m s~ x . Parity selec- 
tion in cyclic dynamos is a subtle issue and space limi- 
tations preclude a thorough discussion here. We only re- 
mark that mean-field BL models have demonstrated that 
the relatively homogeneous nature of tur bulent pumping 
(|Guerrero fc de Gouveia Dal Pind 120081) . the antisym- 
metric structure of a supplementary q-effect operating 
near the base of the convection zone (Dikpati fc Gilmanl 
120011 ). and efficient tu rbulent mixing that peak s in the 
upper convection zone (Ho tta fc Yokovamdl2lH0T ) can all 
help promote negative (dipolar) parity. These mean-field 
results are loosely consistent with the 3D convection sim- 
ulations reported here but warrant a more careful inves- 
tigation which we reserve for future work. 

In summary, we have shown that flux emergence can 
promote cyclic magnetic activity in a convective dy- 
namo by means of the Babcock-Leighton mechanism. 
Although this result is to some extent anticipated by 
mean-field dynamo models, we have demonstrated it for 
the first time in a global convection simulation. The BL 
mechanism is not required to achieve cycles in convective 
dynamo simulations but, as we have shown here, it may 
help shape cycle properties such as the period, ampli- 
tude, and equatorward propagation of toroidal flux. 

However, achieving cyclic activity through the BL 
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mechanism is not easy. To make a significant impact 
on the dynamo, the amplitude of the BL a-effect must 
be comparable to the convective velocity, V c . In partic- 
ular, we find that o.q > 10 m s" 1 is required to induce 
cyclic activity when V c ~ 100 m s _1 in the upper con- 
vection zone. By comparison, typical values of a used in 
mean-field BL dy namo models of the solar cycle are less 
than 1 m s _1 (e.g. Dik pati fc Gilman1l2009tlCnarbonneaul 
l2010h . 

The relatively large value of ao used here, together 
with the strong shear | Vf2| and the efficiency of convec- 
tive transport, can likely account for the very short cycle 
period in this young, rapidly- rotating Sun (f2 = 351©). 
Might other parameters produce a 22-yr period compara- 
ble to the solar cycle? This remains to be seen. Whether 
the cycle period is regulated by a dynamo wave or flux 
transport, the large value of ao needed to promote cyclic 
activity (> 10 m s _1 ) and the short convective turnover 
time scale (~ 20 days) may favor relatively short cycle 
periods. Longer cycle periods can generally be achieved 
in mean-field BL dynamo models b y reducing the effi- 
ciency of poloidal flux tr ansport (e.g lDikpati fc Gilmanl 
12001 ICharbonneaul[20To| ) but we do not have that frce- 
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